Software

In R, there are a number of packages that deal with differential equations. CRAN View

  • deSolve - the main integrator
  • ReacTran - simulating reaction diffusion equations
  • FMA - “flexible modeling
  • bvpSolve - boundary value problems
  • simecol - interactive enviromnent for implementing models

Stochastic Differential Equations have a different set

  • sde - simulatio n and inference of stochastic differential equations
  • GillespieSSA
  • adaptivetau - are for quick simulation

deSolve Package

FME Package

Stands for flexible modeling pacakge, allows one to define some cost function

rootsolver

ODE

Using ode from package deSolve

Example: Prey and Predator

LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator #  \beta x y
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K) # \alpha x
    MortPredator <- rMort * Predator

    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator

    return(list(c(dPrey, dPredator)))
  })
}

pars  <- c(rIng   = 0.2,    # /day, rate of ingestion
           rGrow  = 1.0,    # /day, growth rate of prey
           rMort  = 0.2 ,   # /day, mortality rate of predator
           assEff = 0.5,    # -, assimilation efficiency
           K      = 10)     # mmol/m3, carrying capacity

yini  <- c(Prey = 1, Predator = 2) # initial states of the system
times <- seq(0, 200, by = 1) # initial time
out   <- ode(yini, times, LVmod, pars)
summary(out)
##                Prey    Predator
## Min.      1.0000000   1.8632829
## 1st Qu.   1.9999571   3.9999115
## Median    2.0000000   4.0000000
## Mean      2.0317905   3.9606228
## 3rd Qu.   2.0000751   4.0000418
## Max.      4.2001812   4.9787222
## N       201.0000000 201.0000000
## sd        0.3138139   0.3489079
## Default plot method
plot(out)

## User specified plotting
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "time", ylab = "Conc",
        main = "Lotka-Volterra", lwd = 2)
legend("topright", c("prey", "predator"), col = 1:2, lty = 1:2)

Example: Substrate Producer/Consumer

SPCmod <- function(t, x, parms, input)  {
  with(as.list(c(parms, x)), {
    import <- input(t)
    dS <- import - b*S*P + g*C    # substrate
    dP <- c*S*P  - d*C*P          # producer
    dC <- e*P*C  - f*C            # consumer
    res <- c(dS, dP, dC)
    list(res)
  })
}

## The parameters 
parms <- c(b = 0.001, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)

## vector of timesteps
times <- seq(0, 200, length = 101)

## external signal with rectangle impulse
signal <- data.frame(times = times,
                     import = rep(0, length(times)))

signal$import[signal$times >= 10 & signal$times <= 11] <- 0.2

sigimp <- approxfun(signal$times, signal$import, rule = 2)

## Start values for steady state
xstart <- c(S = 1, P = 1, C = 1)

## Solve model
out <- ode(y = xstart, times = times,
           func = SPCmod, parms = parms, input = sigimp)

## Default plot method
plot(out)

## User specified plotting
mf <- par(mfrow = c(1, 2))

matplot(out[,1], out[,2:4], type = "l", xlab = "time", ylab = "state")
legend("topright", col = 1:3, lty = 1:3, legend = c("S", "P", "C"))
plot(out[,"P"], out[,"C"], type = "l", lwd = 2, xlab = "producer",
  ylab = "consumer")

par(mfrow = mf)

Example: SIR Modeling

A simple example is provided by https://www.r-bloggers.com/2017/01/sir-model-with-desolve-ggplot2/

# definition of the differential equation
sir_diffeq <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dSusceptible <- -beta * Infected * Susceptible
    dInfected <- beta * Infected * Susceptible - gamma * Infected
    dRecovered <- gamma * Infected
    return(list(c(dSusceptible, dInfected, dRecovered)))
  })
}

# solver for differential equation
sir_ode <- function(beta = 1, gamma = .2,
                initial_state = c(Susceptible = 1-1e-6, # initial state defined
                                  Infected = 1e-6,
                                  Recovered = 0),
                times = seq(0, 100, by = 1)) {
  pars <- c(beta = beta, gamma = gamma)
  
  # check that initial state has the correctly named variables
  ode(initial_state, times, sir_diffeq, pars)
}
# Parameter  values
beta <- 1
gamma <- .2

out <- sir_ode(beta = beta, gamma = gamma)

out %>% as.data.frame() %>%
  pivot_longer(Susceptible:Recovered, names_to = "state", values_to = "prop") %>% 
  ggplot(aes(time, prop)) + 
  geom_line(aes(color = state)) +
  theme_minimal() + 
  labs(title = sprintf("SIR, beta = %g, gamma = %g", beta, gamma))

Example: Rossler attractor

This example is a non-linear ordinary differential equation system with chaotic dynamics

\[ \begin{aligned} \frac{dx}{dt} &= -y - z \\ \frac{dy}{dt} &= x + ay \\ \frac{dz}{dt} &= b + z(x -c) \\ \end{aligned} \]

equations x and y are linear which gives the attractor some regularity

rosslerode <- function(t, state, parms) {
  with(as.list(state), {
    dx <- -y - z
    dy <- x + .2 * y
    dz <- .2 + z * (x - 5)
    return(list(c(dx, dy, dz)))
  })
}

yini <- c(x=1, y = 0, z = .9)
times <- seq(0, 100, .05)

out <- ode(times = times, y = yini, func = rosslerode, parms = NULL)
scatterplot3d::scatterplot3d(out[,"x"], out[,"y"], out[,"z"], pch = 19, cex.symbols = .5)

rossler_plotly <- plot_ly(out %>% as.data.frame(), x = ~x, y=~y, z= ~z,
                          marker = list(color= ~time, colorscale = c("#FFFFFF", "#683531")))
rossler_plotly
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Example: Van der Pol Equation

This is a good example of a system that can be either stiff or non-stiff.

  • For large values of \(\mu\), it’s stiff.
  • For small values of \(\mu\), it’s non stiff

The 2nd order ODE can be transformed into a 2 first order equations:

\[ \begin{aligned} y'' - \mu(1 - y^2)y' + y &= 0 \\\\ x' &= y\\ y' &= \mu(1 - x^2)y - x \end{aligned} \]

vdpol <- function(t, state, mu) {
  with(as.list(state), {
    list(c(y, mu * (1 - x^2) *y - x))
  })
}

yini <- c(x = 2, y = 0)
vdpol_stiff <- ode(yini, times = 0:3000, func = vdpol, parms = 1000)
vdpol_nonstiff <- ode(yini, times = seq(0, 30, .01), func = vdpol, parms = 1)

plot(vdpol_stiff, type = "l", main = "IVP ODE, stiff", xlim = c(0, 3000), which = "x")

plot( vdpol_nonstiff,type = "l", main = "IVP ODE, nonstiff", xlim = c(0, 30), which = "x")